#ifndef PROB_RIEMANN_HPP
#define PROB_RIEMANN_HPP
#include "Likelihood.hpp"

class Prob_riemann : public Likelihoods::Probit3 
{
	friend class Likelihoods::Probit3;
	friend class Density::GeomBridge;
	public:
	Prob_riemann(double *Y, mat X, Distribution::Distribution *S, mat s) : Likelihoods::Probit3(Y,X,S,s){
		d=Get_d();
	}
	mat Gt(mat theta)
	{
		mat X=Get_X();
		double *Y=Get_Y();
		int p=X.n_cols;
		mat sum(p,p);
		sum.fill(0);
		double phi_n=Get_phi_n();
		int n=X.n_rows;
		double xb=0;
		Distribution::Distribution *Start=Get_Start();
		for(int i=0;i<n;i++)
		{
			xb=dot(X(i,span::all),theta);
			double gpxb=Phi(xb);
			double pxb=pdf(*d,xb);
			double pxb2=pow(pxb,2);
			double igpxb=gpxb*(1-gpxb);
			
			if(igpxb==0){igpxb=0.00000000001;}
			if(pxb==0){pxb=0.0000000001;}
			if(pxb2==0){pxb2=0.0000000001;}
			//cout << "foo1 "<<foo1;
			//cout << "foo2 "<<foo2;

			mat temp=pxb2/(igpxb)*(X(i,span::all).t()*X(i,span::all));
			sum+=temp;						
		}
		sum=phi_n*(sum-1/((Get_s()(0,0))))+(1-phi_n)*(inv(Start->Get_s()));
	/*	for(int i=0;i<p;i++)
		{
			for(int j=0;j<p;j++)
			{
				if(sum(i,j)<-90000000){
					sum(i,j)=-90000000;
				}else if(sum(i,j)>90000000){
					sum(i,j)=90000000;
				}
				if(sum(i,j)!=sum(i,j))
				{
					cout << "2";
					cout << dot(X(i,span::all),theta) << "\n";
				}
			}
		}
		*/
		
		return sum;

	}
	mat dG(mat theta, int j)
	{
		mat X=Get_X();
		double *Y=Get_Y();
		int p=X.n_cols;
		mat sum(p,p);
		sum.fill(0);
		int n=X.n_rows;
		double phi_n=Get_phi_n();
		double xb=0;
		for(int i=0;i<n;i++)
		{
			xb=dot(X(i,span::all),theta);
			double gpxb=Phi(xb);
			double pxb=pdf(*d,xb);
			double pxb2=pow(pxb2,2);
			double dpxb=-xb*pxb;
			double igpxb2=pow(gpxb*(1-gpxb),2);
			if(igpxb2==0){igpxb2=0.00000001;}
			if(pxb==0){pxb=0.00000001;}
			if(pxb2==0){pxb2=0.00000001;}
			if(dpxb==0){dpxb=0.00000001;}
			//cout << "foo1 "<<foo1;
			//cout << "foo2 "<<foo2;
			//if(pow(gpxb*(1-gpxb),2)==0){cout <<"gpxb=0";}
			//if(pxb!=pxb){cout <<"pxb";}
			//if(pxb2!=pxb2){cout <<"pxb2";}
			double bar0=X(i,j)*pxb*(2*dpxb*gpxb*(1-gpxb));
			double bar2=2*gpxb*pxb2-pxb2;
			double bar1=bar0+bar2;
			double foo1=bar1/igpxb2;
			mat temp=foo1*(X(i,span::all).t()*X(i,span::all));
			sum+=temp;						
		}
		sum=phi_n*sum;
		/*for(int i=0;i<p;i++)
		{
			for(int j=0;j<p;j++)
			{
				if(sum(i,j)<-10000){
					sum(i,j)=-10000;
				}else if(sum(i,j)>10000){
					sum(i,j)=10000;
				}
				if(sum(i,j)!=sum(i,j))
				{
				//	cout << "3";
				//	cout << dot(X(i,span::all),theta) << "\n";
					
				}
			}
		}*/	
	
		return sum;
	}
	private:
	boost::math::normal_distribution<> *d; 


};
#endif
